Non-equilibrium dynamics of Bosonic Mott insulators in an electric field 
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We study the non-equilibrium dynamics of one-dimensional Mott insulating bosons in the presence 
of a tunable effective electric field 8 which takes the system across a quantum critical point (QCP) 
separating a disordered and a translation symmetry broken ordered phase. We provide an exact 
numerical computation of the residual energy Q, the log- fidelity F, the excess defect density D, and 
the order parameter correlation function for a linear-in-time variation of 8 with a rate v. We discuss 
the temporal and spatial variation of these quantities for a range of v and for finite system sizes as 
relevant to realistic experimental setups [J. Simon et al, Nature 472, 307 (2011)]. We show that in 
finite-sized systems Q, F, and D obey Kibble-Zurek scaling, and suggest further experiments within 
this setup to test our theory. 

PACS numbers: 75.10.Jm, 73.43.Nq 



Ultracold atoms in optical lattices provide us with a 
unique opportunity to study both equilibrium phases 
and non-equilibrium quantum dynamics of strongly cou- 
pled bosonic systems near a quantum phase transition 
(QPT) [1]. One system, which has been the subject of 
a recent experimental study, consists of one-dimensional 
(ID) Mott insulating (MI) bosons in the presence of an 
effective electric field £ [2]. It has been shown that this 
system can be described in terms of an effective quantum 
dipole model or, equivalently, an Ising spin model with 
infinitely strong nearest neighbor coupling in the pres- 
ence of both a transverse and a longitudinal field [3] . It 
has also been demonstrated in Ref. [3] , that tuning f to a 
critical value £c leads to a QPT. In the dipole language, 
this transition consists of a change in the ground state 
from a dipole vacuum (which corresponds to n bosons 
at each site) to one with maximal dipoles (which corre- 
sponds to alternating n—1 and n + 1 bosons per site). In 
the spin language, this transition is from the paramagnet 
(PM) to an Ising antiferromagnet (AFM). The interme- 
diate QCP belongs to the Ising universality class [3]. The 
appearance of this AFM order has recently been observed 
using a quantum gas microscope [2] . Theoretical studies 
of the phases of the bosonic Mott insulator in an electric 
field have also been extended to several 2D lattices [4]. 
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FIG. 1: (color online) A pictorial representation of the dipole 
(AFM) and the vacuum (PM) ground states for MI bosons 
with n = 1 across the quantum phase transition. The transi- 
tion occurs at U — 8c — —1.31 (for J = 1). 



The study of non-equilibrium dynamics of closed quan- 
tum systems have seen tremendous progress in recent 
years [5]. One reason for this intense effort has been 
the possibility of experimentally realizing these dynamics 
using ultracold atoms. Indeed, experiments probing non- 
equilibrium phenomena with strongly coupled 2D bosonic 
atoms, well described by the 2D Bose-Hubbard model, 
have been carried out recently [6]. The corresponding 
theoretical studies show a reasonable qualitative match 
with experiments [7]. For the case of ID bosonic MI in 
the presence of an electric field, order parameter dynam- 
ics following sudden quenches of the electric field have 
also been studied theoretically [8]. However, the case of 
finite velocity and, in particular, nearly adiabatic ramps 
of 8 has not been previously explored. 

In this letter, we probe the non-equilibrium dynam- 
ics of the bosonic Mott insulators in the presence of a 
linear-in-time varying electric field (chemical potential 
gradient) 



8(t) =8i^ {£f - 8i)t/r, 



(1) 



where £i and £f are the initial and final values of the 
electric field, and r is the ramp time. We define the 
ramp rate to be v = dt£{t). We look at ramps that 
start from the PM phase with unit boson occupation 
per site and end either in the AFM phase across the 
QPT or at the QCP. Based on the resonant manifold 
model of Ref. [3], we provide an exact numerical com- 
putation for finite system sizes (L), using Exact Diag- 
onalization (ED, L < 16) and time dependent Matrix 
Product States (tMPS, L < 48), of the residual energy 
Q, the log- fidelity F, the number of defects D (a defect, 
for the ramps we address, is a doubly occupied or empty 
lattice site), and the defect correlation function Cij{t). 
We note that the experimental setup of Ref. [2] consti- 
tutes systems of L < 16 lattice sites and measures D as 
well as Cij{t); hence our theoretical analysis constitutes 
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a quantitatively exact description of the dynamics of the 
experimental system providing a valuable guideline to fu- 
ture experiments. Further, for finite sized systems, our 
analysis reveals the manifestation of universal Kibble- 
Zurek-like (KZ) scaling [9-11]. In the past, KZ scaling 
has been shown to work well in integrable systems [12]; 
however, its manifestation in non-integrable systems has 
not been consistently demonstrated [13]. Our work thus 
provides the first realization of KZ scaling in finite-sized 
non-integrable systems which can be tested with an ex- 
isting experimental setup. Furthermore, we investigate 
the crossover from Landau-Zener (LZ) scaling [10, 14] to 
KZ scaling in finite-size systems, and show that it col- 
lapses onto a universal functional form, as suggested in 
Ref. [15]. 

Model — Bosonic atoms in a tilted optical lattice (i.e. 
in an effective time dependent electric field) are well de- 
scribed using the Bose-Hubbard Hamiltonian [3] 

Hb {t) = -to ^ blbj + ^ ^ (n, - 1) - S{t) ^ z • 

{ij) i i 

where (ij) denotes the sum over nearest-neighbors, bi the 
boson annihilation operator at site i, and rii = b\bi the 
boson density operator. Under the condition U^S{t) 
\U — 5(t)|,to, the low-energy dynamics of Hb are cap- 
tured by an effective dipole model [3] 

Hd{t) = [U -£{t)]Y,d\d^-jY,{d\^d^), (2) 

where d£ = bib^j /\/^o(^o + 1) denotes the dipole annihi- 
lation operator which lives on a link i between the neigh- 
boring sites i and j as schematically shown in Fig. 1, 
no the boson occupation of parent Mott insulator, and 
J = to\/^o(^o + 1) the amplitude for creation or an- 
nihilation of a dipole. Henceforth, we shall use the 
units in which = 1, J = 1, and restrict ourselves 
to no = 1. The dipoles satisfy two constraints; first, 
there can be only one dipole on any link £ which renders 
dld£ < 1, and second, two consecutive links can not be 
simultaneously occupied by dipoles dldl_^-^d£-^id£ = 0. 
These constraints render Hd non-integrable; however, 
they also lead to a significant reduction in the Hilbert 
space of Hd which makes ED and tMPS the methods 
of choice for studying the dynamics of the model. In 
addition, we restrict ourselves to studying with peri- 
odic boundary conditions so as to approach the thermo- 
dynamic limit with smaller systems. We note that the 
dipole model can be represented in terms of an Ising- 
like spin-model via the transformation: /S| = 1/2 — djd^, 

Sfy^ = {-i){d£ + {-)dl)/2 [2, 3]. Note that is the 
order parameter for the transition from the PM (dipole 
vacuum) to the AFM (maximal dipole density) state. 

To study the dynamics within exact diagonalization, 
we evolve the wave function using the time dependent 



Hamiltonian of Eq. (2), in which the electric field is tuned 
linearly in time according to Eq. (1) 

ihdtim) = Hd{t)m))- (3) 

We supplement the time evolution with the initial con- 
dition |'0(t = 0)) = IV^g)*, where IV^g)* is the ground 
state wave function of the initial Hamiltonian. Integrat- 
ing Eq. (3) from t = to t = r we obtain the wave 
function at the end of the ramp |?/^(r)). 

We supplement our exact diagonalization studies 
by tMPS, which allows us to study larger system 
sizes. tMPS represents the wave function as |?/^) = 
E{a.}^r^2'''-^n^i'^2...), where are a set of 
X by X matrices indexed by site i and spin (Ti [16]. We 
take advantage of the fact that in the reduced Hilbert 
space the Hamiltonian Hdit)^ Eq. (2), is a sum of single 
site Hamiltonians to perform time evolution. We evolve 
in time via |?/^(t + e)) = P ex.-p{—ieHd{t)) by first exactly 
applying the single site hamiltonian Hd{t) and then pro- 
jecting out configurations with neighboring dipoles using 
the projection operator P. This projection increases the 
MPS bond dimension which is then reduced back to its 
original value x introducing a small error. The method 
becomes exact in the limit e ^ and matrix size x ^ co, 
and its convergence has been numerically checked by ex- 
trapolating in these parameters [17]. Ground states are 
found by evolving the same Hamiltonian to large imagi- 
nary time at fixed S. 

For this work, the specific observables of interest will 
be the residual energy Q, the log-fidelity F, the defect 
number rid (i-e. the number of sites with even parity of 
the boson occupation number), excess defect number 
and the spin-spin correlation function Cij which, at any 
instant are given by 



Q{t) 


= {m\H{t)\m)-EG{t), 


(4) 


Fit) 


= iogmtMGm% 
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nd{t) 


= (^(t)|^^(l + 25|)|^W), 


(6) 


D{t) 


= nrf(t)-{^G(t)lE,(l + 25|)|V'GW), 


(7) 




= {m\s!s^\m), 


(8) 



where Ecit) [|'0g(^))] corresponds to the ground state 
energy [wave function] of the Hamiltonian Hdif). For 
notational brevity, we shall drop the index t from ob- 
servables when evaluating them at the end of the ramp 

(t = T). 

We begin with a discussion of Q and F for finite size 
systems (L < 48) undergoing ramps from the PM phase 
{U-8i = 100) to the QCP {U-£f = U-£c = -1-31) in 
time r. For a finite size system which is always gapped, Q 
and F behave differently than their counterparts in the 
thermodynamic limit (L ^ oo) for which one expects 
KZ-like scaling to manifest in both Q - y{d+z)v/{zv+i) 
and F ~ r^jdv/{zv+i) small v (large r). Here d = 1 is 
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FIG. 2: (color online) Residual energy Q (a) and log- fidelity F (c) as a function of the ramp rate v = dtS{t) for various system 
sizes from 10 < L < 16 (ED) and 24 < L < 48 (tMPS, error bars correspond to time step extrapolation error, bond dimension 
errors are below 10"^). The ramps start in the PM (U — Si — 100) and end at the QCP {U — Sf — —1.31). Dashed lines show 
the indicated power laws. Note the extension of the intermediate Kibble- Zurek- like scaling regime to lower values of v for larger 
system sizes, (b) and (d) show finite size scaling collapse for Q and F for several L. (e) and (f) show Q and F as a function of 
V for ramps starting in the PM {U - Si = 100) and ending in the AFM phase {U - Sf = -100) with 10 < L < 14. Note the 
change in Q power law from v'^ to v'^^'^ . 



the dimensionality, z = 1 is the dynamical critical expo- 
nent, and u = 1 the correlation length critical exponent. 
For very fast ramps, the wave function does not have time 
to evolve during the dynamics, and therefore the behav- 
ior of finite- and infinite-sized systems is similar but not 
universal. As the ramp rate becomes slower {v ~ 1), 
KZ scaling sets in for both infinite- and finite-sized sys- 
tems. However, while KZ scaling is expected to persist 
to infinitely slow ramps in the thermodynamic limit, for 
finite-sized systems it is cut off for ramps slower than a 
critical ramp rate Vc{L) ~ For v < Vc{L)^ Q 

and F are expected to scale as v'^ as in gapped systems 
[10, 14]. These expectations may be formalized in the 
form of scaling functions 

Q ^ L^^(^+^)^/(^^+i)^/^(^L^/^+^), (9) 

F - LV^/^^^+^Vr(^i^^/'^+^), (10) 

where gr{x <C 1) - x2-id+z)u/izu+i) ^ <C 1) - 

^2-du/{zu-\-i) (^ygj.y glow, i.e. LZ regime) and gr{x ^ 
1) ^ const., fr{x ^ 1) ~ const. (KZ regime) [15]. 

The above-mentioned expectations are corroborated in 
panels (a) and (c) of Fig. 2, which show the behavior of Q 
and F as a function of the quench rate for several sys- 
tem sizes 10 < L < 48 on a log-log scale. For both Q and 
F we find scaling for very slow ramps and KZ-like scal- 
ing Q ^ and F ~ v^/^ for intermediate ramps. Finally, 
for fast ramps, we find non power-law behavior followed 



by a plateau. We note that KZ-like scaling only occurs 
in systems of sufficient size (L > 10 for Q and L > 16 for 
F) [18]. The LZ-KZ scaling crossover (Eqs. (9) and (10)) 
is further corroborated in panels (b) and (d) of Fig. 2. 
These plots demonstrate both the scaling collapse of Q 
and F for slower ramps (points to the left) and their devi- 
ations from scaling for faster ramps (points to the right) 
in the non-universal regime. The dashed lines indicate 
the expected form of the scaling functions gr and fr in 
both the LZ and the KZ regimes. Note that for the LZ 
regime, we see collapse even for relatively small system 
sizes (L > 10); however, for the KZ regime collapse only 
occurs for larger system sizes (L > 16 for Q and L > 48 
for F). 

Next, we study ramps that cross the QCP. As the sys- 
tem progresses towards the QCP, the gap closes, and ex- 
citations are produced. At some later time, after passing 
the QCP, the gap reopens, and no new excitations can 
be created. However, although the number of excita- 
tions is conserved in the final adiabatic part of the ramp, 
the energies of individual excitations continue changing. 
Therefore, in ramps across the QCP, we can expect to 
see KZ-like scaling for F but not for Q. These expecta- 
tions are corroborated in panels (e) and (f) of Fig. 2, in 
which we plot Q and F as a function of the ramp rate 
for ramps from PM {U - £i = 100) to AFM phase 
{U-£f = -100). We find KZ-like scaling in F - v^/^ as 



4 




FIG. 3: (color online) (a) Ud/L as a function of the end point of the ramp U — 8f for several ramp rates (all ramps start in 
the PM phase U — 8i = 50). (b) D as a function of v for several L showing KZ scaling (U — Si = 100, U — 8f = —1.31). (c) 
Spin-spin correlation function Ce = {Sf Sij^i) at the end of the ramp as a function I and the ramp rate v for L = 16. 



expected; interestingly Q also scales as v^/^. This hap- 
pens due to the fact that during the adiabatic part of 
the ramp, almost all excitations are converted into singly 
occupied sites, with energy per excitation of ~ \U — £f\. 
Therefore Q becomes identical to the number of excita- 
tion and hence scales in the same way as F. Thus our 
analysis suggest anomalous scaling of Q in ramps across 
the QCP should be interpreted with caution. 

Having obtained the scaling behavior, we concentrate 
on its manifestation for experimentally observable quan- 
tities. We note that the existing experimental setup [6] 
focuses on imaging the parity of the number of bosons on 
each site, after projecting the wave function into a Fock 
state [1, 2, 6]. Effectively, this imaging counts the num- 
ber of empty and doubly occupied sites (i.e. "defects" on 
the PM side). In addition to measuring this defect num- 
ber n^i, the existing experiments can measure the defect- 
defect correlation function gd{^) = + ^)) ^ 
{SfS^j^^) = Ci^i-^i at the end of the ramp. In Fig. 3a, we 
plot the defect density Ud as a function of the U — Sp- 
We note that the final saturation value of rid in the AFM 
phase is a decreasing function of the ramp rate and lies 
between those for a nearly adiabatic ramp and a sudden 
quench. In Fig. 3b, we show that the excess defect num- 
ber D demonstrates similar finite-size KZ scaling as F. 
Since Ud is experimentally measurable, our work demon- 
strates that finite size scaling can be observed within an 
existing experimental setup. Finally, in Fig. 3c, we show 
the behavior of Ci^iJ^i as a function of v and I. We note 
that for slow ramps to the QCP one finds an oscilla- 
tory behavior of indicating the precursor of the 
AFM order present in the critical ground state. As v 
is increased, the system ceases to reach the final ground 
state and these correlations decay. Finally, the state of 
the system at the end of fast ramps is essentially the PM 
ground state, which has no AFM correlations, leading to 
a fiat Ci^iJ^^. These features can be directly picked up in 
experiments, which can serve as a test of our theory. 

In conclusion, we have investigated universal scaling 
dynamics of finite size non-integrable bosonic system fol- 
lowing a finite rate ramp of the effective electric field. 
Our investigation demonstrates two scaling regimes (LZ- 



and KZ-like scaling) with conventional exponents and 
thus differ from prior studies of other non-integrable sys- 
tems [9, 12, 13] which found various anomalous scaling 
exponents. Furthermore, comparing ramps that cross the 
QCP (which show anomalous exponents) to those that 
end at the QCP (which show expected exponents), we 
suggest a possible origin of these anomalous exponents: 
the adiabatic dynamics of the excitations following the 
passage through the quantum critical regime. Finally, 
we compute experimentally measurable quantities such 
as Ud and Cij and demonstrate that the scaling behavior 
studied in this work can be observed in realistic experi- 
ments via measurement of Ud- 
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